On Type Consistency

Currently, the interface contains type parameters for the state and observation objects; while this is ultimately good for forward simulation, we no longer have easy access to the element types. This is used mainly for preallocating particle weights, log evidence, etc and helps avoid unnecessary type conversions. At its worst, this raises errors when facing impossible type conversions.

Taking this to an extreme, I also think the user should keep the model element types consistent throughout the dynamics and observation process. This would allow us to employ a type structure like StateSpaceModel{T, LatentDynamics{Vector{T}}, ObservationProcess{Vector{T}}}, where Base.eltype(::StateSpaceModel) would return T. I would like to reiterate that keeping T consistent is the only way to avoid unnecessary/redundant type promotion.

Awesome stuff. I'll take a look through and make some comments.

Code style is Blue. I use the formatter baked into the Julia VSCode extension set up to format my code on every save. Works like a charm.

A few random thoughts:

• Should resampling take place in its own function? E.g. a step consists of resample, predict, update. For the tracking-in-clutter applications, it seems likely we're going to have to change step to a predict, associate, update loop, so doesn't feel as bad to add extra steps to the particle filter
• On the question of particle storage and ancestry (for naive smoothing), this may be too abstract but I wonder whether it is best left up to the storage object to determine how to save and store this. If the particle filter provides the new states, log weight increments and parent indexes, it can then be down to the storage object to store that as a linked list.
• As mentioned in my message about RTS, I wonder if rather than having a callback function, a callback struct with two methods, post_predict_store and post_update_store would be better. This would allow you to constantly update state in place rather than separating out into proposed_states and filtered_states.

I like the look of the "distribution-focused" version of the Kalman Filter. It might be worth leaving a comment clarifying what that is about for Frederic/Hong since I think we spoke more about that by ourselves. I'm apprehensive about using the name particles. state is much clear to me, where your state can be e.g. a Gaussian, a categorical distribution (HHM), or a weighted collection of point masses (particle collection).

Particle Ancestry

I implemented the particle storage algorithm presented in (Murray, 2015). I include a demonstration at the bottom of this file which informally benchmarks the performance of the algorithm.

@THargreaves since you brought my attention to this algorithm, your comments would be much appreciated. I was aiming for elegance with this commit, and is thus far from optimal in terms of speed. Feel free to scrap anything remotely wasteful.

This looks great and closely matches the initial draft that I made a few months back.

I've pushed a small change that keeps track of the free indices using a stack. On my computer using the parameters in your test script, this reduced the median time from 48ms to 5ms and I assume will have better performance for sequences with fewer dead branches. I haven't tested it in full generality yet though but intuitively this feels it will be worth it most of the time.

Other than that, I think the particle storage code looks great.

I'm not sure about the predict and filter methods dispatching on AncestryContainer though. It would be much more general to have these implementation details abstracted away. Say by a store!(::AbstractParticleContainer, proposed_states idx) or something similar.

I've pushed a small change that keeps track of the free indices using a stack. On my computer using the parameters in your test script, this reduced the median time from 48ms to 5ms...

@THargreaves legendary commit. 10x speed-up in less than 10 lines of changes.

I'm not sure about the predict and filter methods dispatching on AncestryContainer though. It would be much more general to have these implementation details abstracted away. Say by a store!(::AbstractParticleContainer, proposed_states idx) or something similar.

Yeah, it was purely to demonstrate a proof of concept. Now that you fixed the main drawback of the ancestry, I went ahead and added some key-word arguments to let initialise know whether to create an AncestryContainer or a ParticleContianer. It uses a nasty conditional, but I couldn't think of anything else.

I also consolidated my demonstrations to script.jl and moved resampling to a separate file.

A Note on Particle Gibbs

I updated the filtering code such that reference trajectories can be passed as key word arguments. This allows for conditional SMC to be used within particle Gibbs blocks using the highly efficient sparse ancestry storage. My implementation of CSMC is based mostly on Nicolas Chopin's own particles, and is nowhere near as general as the methods defined in AdvancedPS. Again, all suggestions are welcome.

Also relevant to #51, although I'm waiting on the Gibbs extension of AbstractMCMC to consider a serious implementation of Particle Gibbs

I guess one thing of note about moving to kwargs is that we no longer have the fine-grained control of which kwargs are accessible at each time step.

This actually might be a net positive as it allows for a lot more flexibility and avoids the need for extra0 which always felt clunky to me.

Might be worth thinking if there are any scenarios where this could lead to trouble. I can't imagine there is any computational cost involved.

Maybe the closest thing would be that this makes it difficult to use the filter online.

We would write

def simulate(...; dts, kwargs...)
    dt = dts[step]
end

Which then expects us to be appending to a vector of dts as observations come in. Not hard to get around this by creating a "memoryless" vector. E.g.

mutable struct MemorylessVector{T}
    x::T
    i::Int
end

Base.getindex(v:: MemorylessVector, i::Int64) = i == v.i ? x : error(..)